FIELDimageR pipeline:
Intro to R - ‘Drone Imagery Analysis’
Introduction
Agriculture and plant breeding
Images can be used in plant phenotyping to draw inference about many traits:
- Geometric traits (i.e. plant height, leaf area index, lodging, crop canopy cover)
- Canopy spectral texture (spectral features)
- Physiological traits (i.e., chlorophyll, biomass, pigment content, photosynthesis)
- Abiotic/biotic stress indicators (i.e., stomatal conductance, canopy temperature difference, leaf water potential, senescence index)
- Nutrients (nitrogen concentration, protein content)
- Yield
FIELDimageR
FIELDimageR is a R package to analyze images and plant phenotyping, and allows to:
- Crop the image
- Remove soil effect
- Build vegetation indices
- Rotate the image
- Build the plot shapefile
- Extract information for each plot
- Evaluate stand count, canopy percentage, and plant height
FIELDimageR pipeline
1) Example 01
Evaluating multiple FLIGHTS DATES
RGB and DSM - (Potato Breeding) Jeffrey Endelman (UW-Madison)
Required packages
## Installing
# install.packages("sp")
# install.packages("raster")
# install.packages("rgdal")
# install.packages("ggplot2")
# install.packages("agricolae")
# install.packages("reshape2")
# install.packages("devtools")
# install.packages("ggrepel")
# install.packages("lme4")
# install.packages("plyr")
# install.packages("DescTools")
# devtools::install_github("filipematias23/FIELDimageR")
## Necessary packages
library(FIELDimageR)
library(raster)
library(agricolae)
library(reshape2)
library(ggplot2)
library(lme4)
library(plyr)
library(ggrepel)
library(DescTools)
Data (R/FIELDimageR)
- Potato breeding: 140 plots (3x15 feet)
- Population: 97 genotypes
- Six UAV flights: 30 (soil base), 40, 70, and 100 days after planting (DAP)
Basic steps
# Rotating image with theta 2.3:
Test.Rotate<-fieldRotate(Test.Crop, theta = 2.3)
#Test.Rotate<-fieldRotate(Test.Crop, clockwise = F)
Building the plot shapefile
- Dataset - Field notes - phenotype (Download: ‘EX1_Data.csv’)
# Making the field Map
Map<-fieldMap(fieldPlot = Data$Plot,
fieldColumn = Data$Column,
fieldRow = Data$Row,
decreasing = T)
Map - EX1 field design (MAP) should start from bottom to top and left to the right. In other words, plot 1 should be where the blue arrow is showing in the image below.
# Rotating the MAP to correct plots position:
rotate <- function(x) t(apply(x, 2, rev))
Map<-rotate(rotate(Map))
Map# Building the plot shapefile (ncols = 14 and nrows = 10)
# x11()
plotShape<-fieldShape(mosaic = Test.RemSoil$newMosaic,
ncols = 14,
nrows = 10,
fieldData = Data,
ID = "Plot",
fieldMap = Map)# Building indices ("NGRDI","GLI")
Test.Indices<- fieldIndex(mosaic = Test.RemSoil$newMosaic,
Red = 1, Green = 2, Blue = 3,
index = c("NGRDI","GLI"),
myIndex = c("(Red-Blue)/Green"))
Extracting Data
Test.Info<- fieldInfo(mosaic = Test.Indices[[c("NGRDI","GLI","myIndex")]],
fieldShape = plotShape$fieldShape,
buffer = -0.10,
n.core = 3)
Test.Info$fieldShape@data
Estimation of plant height using digital surface model (DSM).
# Uploading files from soil base and vegetative growth:
DSM0 <- stack(paste("./DSM/",DSM[1],sep = ""))
DSM1 <- stack(paste("./DSM/",DSM[3],sep = ""))
plot(DSM1)# Cropping the image using the previous shape:
DSM0.Crop <- fieldCrop(mosaic = DSM0,fieldShape = Test.Crop)
DSM1.Crop <- fieldCrop(mosaic = DSM1,fieldShape = Test.Crop)
# Canopy Height Model (CHM):
DSM0.R <- resample(DSM0.Crop, DSM1.Crop)
CHM <- DSM1.Crop-DSM0.R
plot(CHM)# Rotating the image using the same theta from RGB mosaic:
CHM.Rotate<-fieldRotate(CHM, theta = 2.3)
# Removing the soil using the mask from RGB mosaic:
CHM.RemSoil <- fieldMask(CHM.Rotate, mask = Test.RemSoil$mask)dev.off()
par(mfrow=c(1,3))
# Profile plot:
plot(x = CHM.Draw$Draw1$drawData$x-100230,
y = CHM.Draw$Draw1$drawData$layer,
type="l", col="red",lwd=1,ylim=c(0,1),
xlab="Distance (m)", ylab="EPH (m)")
lines(x = CHM.Draw$Draw2$drawData$x-100230,
y = CHM.Draw$Draw2$drawData$layer,
type="l", col="blue",lwd=1,add=T)
# CHM plot:
plot(CHM.Rotate, col = grey(1:100/100), axes = FALSE, box = FALSE,legend=F)
lines(CHM.Draw$Draw1$drawData$x,CHM.Draw$Draw1$drawData$y, type="l", col="red",lwd=2)
lines(CHM.Draw$Draw2$drawData$x,CHM.Draw$Draw2$drawData$y, type="l", col="blue",lwd=2)
#RGB plot:
plotRGB(Test.Rotate)
lines(CHM.Draw$Draw1$drawData$x,CHM.Draw$Draw1$drawData$y, type="l", col="red",lwd=2)
lines(CHM.Draw$Draw2$drawData$x,CHM.Draw$Draw2$drawData$y, type="l", col="blue",lwd=2)# Extracting the estimate plant height average (EPH):
EPH <- fieldInfo(CHM.RemSoil$newMosaic,
fieldShape = Test.Info$fieldShape,
fun = "mean") #fun="quantile"
EPH$plotValueEvaluating all mosaics in a loop
DataTotal<-NULL
for(i in 2:length(MOSAIC)){
EX1 <- stack(paste("./MOSAIC/",MOSAIC[i],sep = ""))
EX1.Crop <- fieldCrop(mosaic = EX1,
fieldShape = Test.Crop,
plot = F)
EX1.Rotate<-fieldRotate(EX1.Crop,
theta = 2.3,
plot = F)
EX1.RemSoil<-fieldMask(EX1.Rotate, plot = F)
EX1.Indices<- fieldIndex(mosaic = EX1.RemSoil$newMosaic,
Red = 1, Green = 2, Blue = 3,
index = c("NGRDI","BGI", "GLI","VARI"),
myIndex = c("(Red-Blue)/Green"),
plot = F)
EX1.Info<- fieldInfo(mosaic = EX1.Indices[[c("NGRDI","BGI", "GLI","VARI","myIndex")]],
fieldShape = plotShape$fieldShape,
n.core = 3)
DSM0 <- stack(paste("./DSM/",DSM[1],sep = ""))
DSM1 <- stack(paste("./DSM/",DSM[i],sep = ""))
DSM0.Crop <- fieldCrop(mosaic = DSM0,fieldShape = Test.Crop,plot = F)
DSM1.Crop <- fieldCrop(mosaic = DSM1,fieldShape = Test.Crop,plot = F)
DSM0.R <- resample(DSM0.Crop, DSM1.Crop)
CHM <- DSM1.Crop-DSM0.R
CHM.Rotate<-fieldRotate(CHM, theta = 2.3,plot = F)
CHM.RemSoil <- fieldMask(CHM.Rotate, mask = EX1.RemSoil$mask,plot = F)
EPH <- fieldInfo(CHM.RemSoil$newMosaic,
fieldShape = EX1.Info$fieldShape,
fun = "mean")
DataTotal<-rbind(DataTotal,
data.frame(DAP=as.character(do.call(c,strsplit(MOSAIC[i],split = "_"))[2]),
EPH$fieldShape@data))
# Making map plots
fieldPlot(fieldShape=EPH$fieldShape,
fieldAttribute="NGRDI",
mosaic=EX1.Rotate,
color=c("red","green"),
# min.lim = 0,
# max.lim = 0.35,
alpha = 0.5,
round = 2)
print(paste("### Completed: ", "Mosaic_",i," ###",sep=""))
}
# Organizing the data table:
colnames(DataTotal)<-c(colnames(DataTotal)[-c(dim(DataTotal)[2])],"EPH") # layer=EPH
DataTotal<-DataTotal[,!colnames(DataTotal)%in%c("ID","ID.1","PlotName")] # Removing column 12 ("ID.1")
DataTotal
# Saving extracted data "DataTotal.csv":
write.csv(DataTotal,"DataTotal.csv",row.names = F,col.names = T)
- Graphics
DataTotal$Name<-as.factor(as.character(DataTotal$Name))
DataTotal$Row<-as.factor(as.character(DataTotal$Row))
DataTotal$Column<-as.factor(as.character(DataTotal$Column))
DataTotal$DAP<-as.numeric(as.character(DataTotal$DAP))
DataTotal$NGRDI<-as.numeric(as.character(DataTotal$NGRDI))
DataTotal$BGI<-as.numeric(as.character(DataTotal$BGI))
DataTotal$EPH<-as.numeric(as.character(DataTotal$EPH))
ggplot(DataTotal, aes(x = NGRDI,fill=as.factor(DAP))) +
geom_density(alpha=.5,position = 'identity') +
facet_wrap(~DAP,ncol = 1)+
scale_fill_grey(start=1, end=0)+
labs(y="#genotypes",x="NGRDI", fill="DAP") +
theme_bw()
- Linear Regression
DataTotal.Reg<-subset(DataTotal,DAP=="40")
DataTotal.Reg$Check<-as.character(DataTotal.Reg$Name)
DataTotal.Reg$Check[!DataTotal.Reg$Check%in%c("G43","G44","G45")]<-""
ggplot(DataTotal.Reg,aes(y=NGRDI, x=EPH)) +
geom_point() +
geom_smooth(method=lm)+
labs(y="EPH",x="NGRDI",fill="",alpha="")+
geom_vline(aes(xintercept=0.02),col="red", linetype = 2, size=0.7) +
theme_bw()+
geom_text_repel(aes(label = Check),
size = 3.5, col="red",
box.padding = 0.35,
point.padding = 0.5,
segment.color = 'grey50')+
theme(legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
- Heritability
DAP<-unique(DataTotal$DAP)
H2<-NULL
for(h in 1:length(DAP)){
mod<-lmer(NGRDI~Row+Column+(1|Name),DataTotal[as.character(DataTotal$DAP)==DAP[h],])
H2.a<-c("NGRDI",DAP[h],as.data.frame(VarCorr(mod))$vcov[1]/sum(as.data.frame(VarCorr(mod))$vcov))
mod<-lmer(EPH~Row+Column+(1|Name),DataTotal[as.character(DataTotal$DAP)==DAP[h],])
H2.b<-c("EPH",DAP[h],as.data.frame(VarCorr(mod))$vcov[1]/sum(as.data.frame(VarCorr(mod))$vcov))
H2<-rbind.data.frame(H2,rbind(H2.a,H2.b))
}
colnames(H2)<-c("Trait","DAP","H2")
H2$H2<-as.numeric(as.character(H2$H2))
H2$DAP<-as.numeric(as.character(H2$DAP))
ggplot(H2,aes(x=as.factor(DAP),y=H2,fill=as.factor(DAP)))+
geom_bar(stat="identity")+
facet_wrap(~Trait)+
scale_fill_grey(start=0.8, end=0.2)+
labs(x="Days After Planting (DAP)", fill="")+
geom_text(aes(label=round(H2,2)), vjust=1.6, color="white", size=6)+
theme_bw()+
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5))
Area under the curve (AUC)
DataTotal1<-DataTotal[as.character(DataTotal$Name)%in%c("G43","G44","G45"),]
ggplot(data=DataTotal1, aes(x=as.numeric(DAP), y= NGRDI, col= Name, group=Name)) +
geom_point(size=6)+
geom_line(size=1.2) +
scale_color_grey(start=0.8, end=0.2)+
labs(x="Days After Planting (DAP)", fill="", col="")+
theme_linedraw()Trait<-c("NGRDI") # c("GLI","EPH")
Plot<-as.character(unique(DataTotal$Plot))
DataAUC<-NULL
for(a1 in 1:length(Plot)){
D1<-DataTotal[as.character(DataTotal$Plot)==Plot[a1],]
x1<-c(0,as.numeric(D1$DAP))
y1<-c(0,as.numeric(D1[,Trait]))
DataAUC <- rbind(DataAUC,
c(NGRDI_AUC=AUC(x = x1[!is.na(y1)], y = y1[!is.na(y1)]),
Name=unique(as.character(D1$Name)),
Trait=unique(D1$Trait),
Row=unique(D1$Row),
Column=unique(D1$Column)))}
DataAUC<-as.data.frame(DataAUC)
DataAUC$NGRDI_AUC<-as.numeric(as.character(DataAUC$NGRDI_AUC))
DataAUC$Name<-as.factor(DataAUC$Name)
DataAUC$Row<-as.factor(DataAUC$Row)
DataAUC$Column<-as.factor(DataAUC$Column)
DataAUC - AUC Heritability
# Mixed model:
mod<-lmer(NGRDI_AUC~Row+Column+(1|Name),DataAUC)
H2<-as.data.frame(VarCorr(mod))$vcov[1]/sum(as.data.frame(VarCorr(mod))$vcov)
H2
ggplot(DataAUC, aes(x = NGRDI_AUC)) +
geom_histogram(aes(y=..density..), colour="black", fill="white")+
geom_density(alpha=.5,position = 'identity', fill="cadetblue") +
labs(y="#genotypes",x=paste("Area under the curve (AUC_NGRDI: H2=",round(H2,2),")",sep="")) +
theme_bw()
2) Example 02
Evaluating multiple FIELD TRIALS
RGB and DSM - (Potato Breeding) Jeffrey Endelman (UW-Madison)
### Orthomosaic using Open Drone Map (ODM)
Follow the OpenDroneMap’s documentation according to your operating system (Windows, macOS or Linux) to install WebODM.
Start WebODM and +Add Project to upload the RGB images from the Tropical Forage Breeding - Embrapa Beef Cattle:
- Donwload Example 1:
Data_ODM
After the running process ‘completed’, download the odm_orthophoto.tif and dsm.tif to upload in R. Then follow the pipeline of FIELDimageR.
### Data extraction (R/FIELDimageR)
#### Required packages
# Installing
install.packages("sp")
install.packages("raster")
install.packages("rgdal")
install.packages("readxl")
install.packages("ggplot2")
install.packages("agricolae")
install.packages("reshape2")
install.packages("devtools")
install.packages("lme4")
install.packages("plyr")
devtools::install_github("filipematias23/FIELDimageR")
# Necessary packages
library(FIELDimageR)
library(raster)
library(readxl)
library(agricolae)
library(reshape2)
library(ggplot2)
library(lme4)
library(plyr)
Donwload:
odm_orthophoto.tif
Evaluating each trial per time:
- Dataset - Field notes - phenotype (Download: ‘EX2_Data.xlsx’)
# Dataset - Field notes - phenotype (.exel or .txt):
Data1<-Data<-read_excel('EX2_Data.xlsx',1)
# Field map ID identification (ID="Plot"):
Map1<-fieldMap(fieldPlot = Data1$Plot,fieldColumn = Data1$Column, fieldRow = Data1$Row,decreasing = T)
Map1 # Building the plot shapefile (ncols = 14 and nrows = 10)
x11()
EX2.Shape.1<-fieldShape(mosaic = EX2.Crop.1, ncols = 14, nrows = 10, fieldData = Data1, ID = "Plot", fieldMap = Map1)
# Trial: 02
EX2.Crop.2 <- fieldCrop(mosaic = EX2.RemSoil$newMosaic)
Data2<-Data<-read_excel('EX2_Data.xlsx',2)
Map2<-fieldMap(fieldPlot = Data2$Plot,fieldColumn = Data2$Column, fieldRow = Data2$Row,decreasing = T)
EX2.Shape.2<-fieldShape(mosaic = EX2.Crop.2, ncols = 14, nrows = 9, fieldData = Data2, ID = "Plot", fieldMap = Map2) # Trial: 03
EX2.Crop.3 <- fieldCrop(mosaic = EX2.RemSoil$newMosaic)
Data3<-Data<-read_excel('EX2_Data.xlsx',3)
Map3<-fieldMap(fieldPlot = Data3$Plot,fieldColumn = Data3$Column, fieldRow = Data3$Row,decreasing = T)
EX2.Shape.3<-fieldShape(mosaic = EX2.Crop.3, ncols = 14, nrows = 13, fieldData = Data3, ID = "Plot", fieldMap = Map3) # Trial: 04
EX2.Crop.4 <- fieldCrop(mosaic = EX2.RemSoil$newMosaic)
Data4<-Data<-read_excel('EX2_Data.xlsx',4)
Map4<-fieldMap(fieldPlot = Data4$Plot,fieldColumn = Data4$Column, fieldRow = Data4$Row,decreasing = T)
EX2.Shape.4<-fieldShape(mosaic = EX2.Crop.4, ncols = 14, nrows = 10, fieldData = Data4, ID = "Plot", fieldMap = Map4)
# Combining shapefiles
EX2.Shape<-rbind(EX2.Shape.1$fieldShape,EX2.Shape.2$fieldShape,EX2.Shape.3$fieldShape,EX2.Shape.4$fieldShape)
plot(EX2.Shape) # Building indices ("NGRDI","BGI","GLI","SCI")
EX2.Indices<- fieldIndex(mosaic = EX2.RemSoil$newMosaic,
index = c("NGRDI","BGI"),
myIndex = c("(Red-Blue)/Green")) # Extracting Data
EX2.Info<- fieldInfo(mosaic = EX2.Indices[[c("NGRDI","BGI","myIndex")]],
fieldShape = EX2.Shape,
buffer = -0.05,
n.core = 3)
NewData<-EX2.Info$fieldShape@data
NewData # Making map plots
fieldPlot(fieldShape=EX2.Info$fieldShape,
fieldAttribute="NGRDI",
mosaic=EX2.Indices,
color=c("red","blue"),
alpha = 0.4,
round = 2)
Estimation of plant height using a high throughput phenotyping platform. (Donwload:
dsm.tif)# Extracting the estimate plant height average (EPH): EPH <- fieldInfo(DSM.S$newMosaic, fieldShape = EX2.Info$fieldShape, fun = "mean",n.core = 3) EPH$plotValue# Correlation NewData2<-EPH$fieldShape@data[,c("dsm","NGRDI","BGI","myIndex")] cor1<-correlation(NewData2) cor1$correlation cor1$pvalue# Regression NewData3<-EPH$fieldShape@data[,c("Trial","Name","Block","Row","Column","dsm","NGRDI")] NewData3$Trial<-as.factor(NewData3$Trial) NewData3$Name<-as.factor(NewData3$Name) NewData3$Block<-as.factor(NewData3$Block) NewData3$Row<-as.factor(NewData3$Row) NewData3$Column<-as.factor(NewData3$Column) NewData3$dsm<-as.numeric(as.character(NewData3$dsm)) NewData3$NGRDI<-as.numeric(as.character(NewData3$NGRDI)) ggplot(NewData3,aes(x=dsm,y=NGRDI, col=Trial))+ facet_wrap(~Trial, scales = "fixed",ncol = 1)+ geom_point() + geom_smooth(method=lm)+ theme_bw()fieldPlot(fieldShape=EPH$fieldShape, fieldAttribute="dsm", mosaic=EX2.Indices, color=c("red","blue"), alpha = 0.4, round = 2)- Heritability
# NGRDI
str(NewData3)
mod.T1<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T1",])
H2.T1<-as.data.frame(VarCorr(mod.T1))$vcov[1]/sum(as.data.frame(VarCorr(mod.T1))$vcov)
mod.T2<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T2",])
H2.T2<-as.data.frame(VarCorr(mod.T2))$vcov[1]/sum(as.data.frame(VarCorr(mod.T2))$vcov)
mod.T3<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T3",])
H2.T3<-as.data.frame(VarCorr(mod.T3))$vcov[1]/sum(as.data.frame(VarCorr(mod.T3))$vcov)
mod.T4<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T4",])
H2.T4<-as.data.frame(VarCorr(mod.T4))$vcov[1]/sum(as.data.frame(VarCorr(mod.T4))$vcov)
NewData4<-data.frame(Trail=factor(c("T1","T2","T3","T4")),
H2=as.numeric(c(H2.T1,H2.T2,H2.T3,H2.T4)))
ggplot(NewData4,aes(x=Trail,y=H2,fill=Trail))+
geom_bar(stat="identity")+
geom_text(aes(label=round(H2,2)), vjust=1.6, color="white", size=6)+
# Estimated plant height
mod.T1<-lmer(dsm~Row+Column+(1|Name),NewData3[NewData3$Trial=="T1",])
H2.T1<-as.data.frame(VarCorr(mod.T1))$vcov[1]/sum(as.data.frame(VarCorr(mod.T1))$vcov)
mod.T2<-lmer(dsm~Row+Column+(1|Name),NewData3[NewData3$Trial=="T2",])
H2.T2<-as.data.frame(VarCorr(mod.T2))$vcov[1]/sum(as.data.frame(VarCorr(mod.T2))$vcov)
mod.T3<-lmer(dsm~Row+Column+(1|Name),NewData3[NewData3$Trial=="T3",])
H2.T3<-as.data.frame(VarCorr(mod.T3))$vcov[1]/sum(as.data.frame(VarCorr(mod.T3))$vcov)
mod.T4<-lmer(dsm~Row+Column+(1|Name),NewData3[NewData3$Trial=="T4",])
H2.T4<-as.data.frame(VarCorr(mod.T4))$vcov[1]/sum(as.data.frame(VarCorr(mod.T4))$vcov)
NewData4<-data.frame(Trail=factor(c("T1","T2","T3","T4")),
EPH=as.numeric(c(H2.T1,H2.T2,H2.T3,H2.T4)))
ggplot(NewData4,aes(x=Trail,y=EPH,fill=Trail))+
geom_bar(stat="identity")+
geom_text(aes(label=round(EPH,2)), vjust=1.6, color="white", size=6)+
theme_bw()Reference
Matias, F.I. (2019). FIELDimageR Pipeline (tutorial). https://github.com/filipematias23/FIELDimageR
Matias, F.I.; Caraza-Harter, M.; Endelman, J.B. (2020). FIELDimageR: A R Package to Analyze Orthomosaic Images from Agricultural Field Trials. The Plant Phenome Journal. DOI:10.1002/ppj2.20005